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ABSTRACT 

. . . Formation of massive stars by accretion requires a high accretion rate of A/, > 10^^ Mq/yt to 

OO I overcome the radiation pressure barrier of the forming stars. Here, we study evolution of protostars 

. accreting at such high rates, by solving the structure of the central star and the inner accreting envelope 

' simultaneously. The protostellar evolution is followed starting from small initial cores until their arrival 

, at the stage of the Zero-Age Main Sequence (ZAMS) stars. An emphasis is put on evolutionary features 

Ch ' different from those with a low accretion rate of Af, ~ 10~^ Mq/jt, which is presumed in the standard 

■ ^ ■ scenario for low-mass star formation. With the high accretion rate of Af* 10~^ Mq/jt, the protostellar 

radius becomes very large and exceeds 100 Rq. Unlike the cases of low accretion rates, deuterium burning 
hardly affects the evolution, and the protostar remains radiative even after its ignition. It is not until 
the stellar mass reaches ~ 40 Af© that hydrogen burning begins and the protostar reaches the ZAMS 
phase, and this ZAMS arrival mass increases with the accretion rate. These features are similar to those 
f~i , of the first star formation in the universe, where high accretion rates are also expected, rather than to 

Qh' the present-day low-mass star formation. At a very high accretion rate of > 3 x 10"'^ Mq/jt, the total 

JL I luminosity of the protostar becomes so high that the resultant radiation pressure inhibits the growth of 

5^ . the protostars under steady accretion before reaching the ZAMS stage. Therefore, the evolution under 

' the critical accretion rate 3 x 10"'^ Mq/jt gives the upper mass limit of possible pre-main-sequence 

^ I stars at ~ 60 Mq. The upper mass limit of MS stars is also set by the radiation pressure onto the 

dusty envelope under the same accretion rate at ~ 250 A^©. We also propose that the central source 
^ ' enshrouded in the Orion KL/BN nebula has effective temperature and luminosity consistent with our 

^ , model, and is a possible candidate for such protostars growing under the high accretion rate. 

^4 I Subject headings: accretion - stars: early-type - stars: evolution - stars: formation - stars: 

i pre-main-sequence 
^ ■ 

■rj" ' 1. INTRODUCTION 

' Massive (> 8 A/©) stars make significant impacts on the interstellar medium via various feedback processes:, e.g., UV 
radiation, stellar winds, and supernova explosions. Such feedback processes sometimes trigger or regulate nearby star- 
25 formation activity. Their thermal and kinetic effects are important factors in the phase cycle of the interstellar medium. 
Strong coherent feedback by many massive stars causes galactic-scale dynamical phenomena, such as galactic winds. 
Furthermore, massive stars dominate the light f rom distant galaxies. Th e cosmic star formati on rate is mostly estimated 
with the hght observed from massive stars (e.g.. iMadau et al.. 1996. : Hopkins fc Beacomll2006[ ). Despite such importance, 
. , the question of the formation process of massive stars still remains open. As for lower- mass stars, there is a widely-accepted 
. formation scen ario, where gravitation al collapse of molecular cloud cores leads to subsequent mass accretion onto tiny 
protostars fe.g.. lShu. Adams fc Lizaiio 1987). If one directly applies this scenario to massive star forma tion, however, some 
difficulties arise in the main accretion phase (e.g., see lZinnecker fc Yorkel[2007t iMcKee fc Ostrikeil2007l for recent reviews). 
The main difficulty in the formation of massive stars is the very strong radiation pressure acting on a dusty envelope. The 
radiative repulsive effect becomes quite strong at the dust destruction front, where the accretion fiow gets much of the 
outward momentum of radiation. Therefore, a necessary condition for massive star formation is to overcome this barrier 
at the dust destruction front. However, some theoretical work has shown t hat the accretion rate of M ^, ^ 10~^ Mq/yi 
expected of lower- mass protostars is too low to overcome the barrier (e.g., IWolfire fc Cassineliilll987l hereafter WC87). 
Because of this diffi culty, various scen arios different from the standard accreti on paradigm have been p roposed such as 
stellar mergers (e.g.. IStahler. Palla fc Ho 2000) and competitive accretion (^e.g.. iBonnell. Bate fc Zinneck er 1998). 

An origin of this disputed situation is uncertainty concerning the initial condition for massive star formation. This 
uncertainty partly originates from the scarcity of massive stars; most of the massive star forming regions are distant. In 
addition, the initial condition is easily disrupted by feedback from newly formed massive stars. Since the accretion rate 
reflects the thermal state of the original molecular core, it is still uncertain which rates should be favored for growing 
high-mass protostars. In fact, some observations of young high-mass sources suggest high accretion rates. Signatures 
of infall motion are detected with various line observations toward high-mass protostellar objects (HMPOs) and hyper- 
/ultra-compact H II regions, and the derived accretion rates are 10-4 - 10-3 Af©/yr (e.g.. iFuller. Williams fc SridharanI 
l2005t Beltran et al] 20061: Keto fc Woodl[200^ . S imilar high accret i on rates are also inferr ed by SED fitting of hot cores 
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(jOsorio, Li zano fc D'Alessiol 119990 and HMPOs ll Fazal et al.]|2007t iKumar fc Gravel l2008f) . Molecular outflows are also 



ubiquitous in high-mass star-forming regions, and estimated high mass outflow rates also suggest high accretion rates 
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(e.g.. iBeuther et al.l 20021: Zhang et alj 200^. Such hig h accretion rate s have the adv a ntage of overcoming the radiation 

i lKahnlHOTl 



pressure barrier ([Larson k. Starrfieldlll97H iKahnI 11974) . Theoretically, iNakano et al.l (j2000l ) have suggested a protostar 



growing at the very high accretion rate of ^ 10~^ Mp,/ yr to explain the low radiation temperature of scattered light from 
Orion KL region (|Morino et al. 199^). iMcKee fc TanI (2002, 20031 have predicted that massive molecular cloud cores, 
from which a few massive stars form, should be domin ated by supersonic turbulence, env isioning that such cores, if they 
exist, are embedded in a high-pressure environment. iKrumholz. Klein fc McKe3 (j2007[ ) have simulated collapse of the 
turbulent cores performing radiation- hydrodynamical calculations, and demonstrated that the accretion rates attain more 
than 10"'' Mq/ji in their simulations. Some recent observations are getting a close lo ok at the initial condition of massive 
star formation, an d have found strong candidates for high-mass pre-stellar cores fe.g.. lRathborne. Simon fc Jackson, 2C)07l : 
iMotte et al.ll2007l ). Further detailed observations will verify the scenario. 

The main targets of this paper are protostars growing at the high accretion rates of M^, > 10~^ Mq/jt. Detailed 
modeling of such protostars should be useful for future high-resolution observations and simulations of massive star 
formation. However, previous studies on such protostars are fairly limited. The protostellar evolution has been well 
studied for low-mass (< 1 Mq) and intermediate-mass (< 8 Mq) protostars by detailed numerical calculations solving 
the stellar structure, but these studies have focused on evolution at the low accretion rate of ~ 10~^ Mq/jv (e.g.. Stabler, 
Shu fc Taam 1980a,b, hereafter SST80a,b; Palla fc Stabler 1990, 1991, hereafter PS91, 1992; Beech fc Mitalas 1994) 
Maeder and coworkers have calcul ated the protostellar evolution under accretion rates growing with the stellar mass, 
which finally exceeds 10"'' Mq/jt (jBernasconi fc Maeded [19961 : iBehrend fc Maeded[200]] : iNorberg fc Maedeil[2000l) . but 
it is still uncertain how the evolution changes with accretion rates, and which physical mechanisms cause differences. 
Some authors have used polytropic one-zone models originally invented for l ow accretion rates even at such a high rate 
(e.g., INakano et all 120001 : iMcKee fc Tm! 120031 : IKrumholz fc Thompson! 120071 ). but its validity remains obscure owing to 
the lack of detailed cal culations. By coincidence, the similar high accretion rates are expected for protostars forming in 
the early universe (e.g.. lOmukai fc Nishilll998l : lYoshida et al.|[2006D . This simp ly reflects high temperatures in p rimordial 
pre-stellar clou ds. Evolu tion of such primordial protostars has been studied by Istahle r. Palla fc SalpeteJ (jl986l . hereafter 
SPS86) and O mukai fc P alla (2001, 2003). Therefore, our other motivation is to investigate how the protostellar evolution 
changes with metallicities. 

To summarize, our goal in this paper is to answer the following questions; 

• What are the properties of massive protostars (e.g., radius, luminosity, and effective temperature) growing at high 
accretion rates? Can we observe any signatures of their characteristic properties? 

• How different are high-mass protostars from low- and intermediate-mass ones, for which lower accretion rates are 
presumed? How about the difference from the protostellar evolution of primordial stars? Furthermore, if the 
protostellar evolution significantly varies with accretion rates and metallicities, what causes such differences? 

• What are the consequences of protostellar evolution with high accretion rates for the formation and feedback 
processes of massive stars? How massive can a star get before its arrival at the Zero- Age Main Sequence (ZAMS)? 
What is the maximum mass of stars that can be formed? 

In order to answer these questions, we solve the stru cture of accreting protostars with various accretion rates and metal- 
licities with detailed numerical calculations (also see I Yorke^ Bodcnhcimcr 2008, for a similar effort). 

The organization of this paper is as follows: In § |2l we briefly review the basic procedure to construct numerical models 
of accreting protostars and their surrounding envelopes. The subsequent § |3| is the main part of this paper, where our 
numerical results are presented. First, we investigate the protostellar evolution of two fiducial cases with the accretion 
rates of 10~^ and 10~^ Mq /yi in 13.11 and 13.21 After that, we show more general variations of protostellar evolution, 
e.g. , mass accretion rates in § 13.31 and metallicities in § 13.41 In § 13.51 we present the protostellar evolution at very high 
accretion rates exceeding 10"'^ Mq/ji, and show that the steady accretion is limited by a radiation pressure barrier acting 
on a gas envelope. We further discuss some implications based on our numerical results. In §|4l we examine which feedback 
effects can affect the growth of protostars for a given accretion rate. We also explore some observational possibilities of 
detecting signatures of the supposed high accretion rates in §|5] Finally, §111 is assigned to summary and conclusions. 

2. NUMERICAL MODELING 

2.1. Outline 

We employ the method of calculation developed by SPS86 and PS91. In this section, only the outline of modeling is 
overviewed. A detailed description of the method is presented in Appendix IA1 

As shown in Figured! the whole system can be divided into two parts of different nature, i.e., an accreting envelope and 
a hydrostatic core. This hydrostatic core is also called a protostar since it eventually grows into a star by accretion. The 
outer part of the accreting envelope contains dust grains, while the warm (> 2000K) part near the protostar is dust-free 
as a result of their evaporation. Despite the absence of the dust, the innermost dense part of the envelope becomes 
optically thick to gas continuum opacity and the photosphere appears outside the accretion shock front in the case of a 
high accretion rate. Such an inner envelope is called the radiative precursor. We here consider only the protostar and 
the radiative precursor under a given constant accretion rate (see Fig. [T|). The dusty outer envelope is not included in 
our formulation, although stellar feedback exerted on the envelope may be important in the last stage of accretion: the 
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Fig. 1. — A schematic figure of a protostar and surrounding accretion flow. The accretion shock forms at the stellar surface. When the 
flow becomes optically thick before hitting on the shock, the photosphere forms outside the stellar surface. The optically thick part of the 
accreting envelope is called the radiative precursor. A dust cocoon surrounds the protostar at larger radius. The inner boundary of the dust 
cocoon is a dust destruction front. Most of the light from the protostar is absorbed at the dust destruction front once, and reemitted as 
infrared light. In our calculations, we solve the detailed structure only of the protostar and accretion flow far inside the dust destruction 
front, which are enclosed by the thick dashed line. We discuss possible feedback effects on the outer dust cocoon in §|4] 



mass accretion can be terminated finally by radiation pressure or other protostellar feedback processes (e.g., stellar wind) 
exerted on the dusty envelope. For the present, we presume that the protostar continues to grow at a given mass accretion 
rate and look for solutions of protostars with steady-state accretion. Discussion on feedback that halts the accretion is 
deferred to § SI 

We calculate protostellar evolution by constructing a time sequence of quasi-steady structures of the protostar and 
accreting envelope. We solve the stellar structure equations for the protostar. For the accreting envelope, we adopt 
different treatments depending on whether or not it is opaque to the gas opacity. If the flow remains optically thin and no 
radiative precursor exists, the free-fall is assumed to be outside the star. For the opaque flow, on the other hand, we solve 
the structure of the radiative precursor by using the equations for a steady-state flow. At the outer boundary, which is 
taken to be at the photosphere, the flow is assumed to be in free fall. After solving the protostar and accreting envelope 
individually, these solutions arc connected at the accretion shock front by the radiative shock condition (e.g., SSTSOa). 
The shooting method is adopted for solving radial structure. This procedure is repeated until the required boundary 
conditions are satisfied. 

We start calculation from a very small protostar, typically M*_o = O.OIAIq. The initial models are constructed following 
SSTSOb (see Appendix lA.2l for detail). Although this choice is rather arbitrary, the star converges immediately to a certain 
structure appropriate for accretion at a given rate. Specific conditions of the initial models do not affect the evolution 
thereafter. The evolution is followed by increasing the stellar mass owing to accretion step by step. This procedure is 
repeated until the star reaches the ZAMS phase after the onset of hydrogen burning. In a few runs with very high accretion 
rates, however, steady accretion becomes impossible before arrival at the ZAMS and the calculation is terminated at this 
moment. 

2.2. Calculated Runs 

The calculated runs and their input parameters are listed in Table. 1. Cases with a wide range of accretion rates are 
studied, starting from 10~^ Mq/jt up to 6 x 10~^ Mq/yt. The lowest values of 10~^ — 10~^ Mq/jt are typical for 
low-mass star formation, while high rates > lO^'* Mq/jt are those envisaged in the accretion scenario of massive star 
formation. The initial deuterium abundance of [D/H] = 2.5 x 10^^ is adopted as a fiducial value, following previous work 
(e.g., SSTSOa, PS91) for comparison. To assess the role of deuterium burning, we also calculate "noD" runs, where the 
deuterium is absent. Most of the runs arc for the solar metallicity ^0(^0. 02). Two runs with lO^'^ and lO^^MQ/yr are 



4 T. Hosokawa & K. Omukai 



Table 1. Calculated Runs and Input Parameters 



Run 


(Me/yr)- 


[D/H] (10-5)^ 








reference-^ 


MD6x3 


6 X 10-3 


2.5 


0.02 


0.5 


43.9 


§1131 


MD4x3 


4 X 10"3 


2.5 


0.02 


0.3 


33.0 


§[33] 


MD3x3 


3 X 10-3 


2.5 


0.02 


0.2 


27.8 


§[S3| 


MD3x3-zO 


3 X 10-3 


2.5 


0.0 


0.2 


26.4 


§[33| 


MD3 


10-3 


2.5 


0.02 


0.05 


15.5 


§133 ESI, [311 [S3| 


MD3-noD 


10-3 


0.0 


0.02 


0.05 


15.5 


oco 


MD3-Z0 


10-3 


2.5 


0.0 


0.05 


13.0 


§[311 


MD4 


10-4 


2.5 


0.02 


0.01 


7.9 


§[33] 


MD4-noD 


10-4 


0.0 


0.02 


0.01 


7.9 


§[33] 


MD5 


10-5 


2.5 


0.02 


0.01 


3.7 


S [QIO ■ lO 


MD5-noD 


10-5 


0.0 


0.02 


0.01 


3.7 


§[321 [33] 


MDS-dhl 


10-5 


1.0 


0.02 


0.01 


3.7 


Appendix IB. 11 


MD5-dh3 


10-5 


3.0 


0.02 


0.01 


3.7 


Appendix IB. 11 


MD5-Z0 


10-5 


2.5 


0.0 


0.01 


2.9 


§[311 


MD5-ps919 


10-5 


2.5 


0.02 


1.0 


4.2 


Appendix IB. 21 


MD6 


lo-'^ 


2.5 


0.02 


0.01 


1.6 


§[33] 


MD6-noD 


10-6 


0.0 


0.02 


0.01 


1.6 


§[33] 



a : mass accretion rate, h : initial number fractional abundance of deuterium, c : metallicity, d : mass of initial core 
model, e : radius of initial core model, / : subsections where numerical results of each run are presented, g : initial model 
is the same as PS91 



calculated for the metal-free gas (runs MD3-zO and MD5-zO) to see the effects of different metallicities. For 10-5 Mq/jt, 
variations in deuterium abundances (runs MD5-dhl and MD5-dh3) and in initial models (MD5-ps91) are studied and 
presented in Appendix [B] for comparison with previous calculations. 

Initial stellar mass in each run is taken to be sufhciently small, typically A/»^o = 0.01 M©. For high accretion rates 
M, > 10-3 Mo/yr, somewhat more massive initial stars (see Table. 1) are used since convergence of calculation was not 
achieved for lower mass ones. The radii of initial models are also listed in Table. 1. 

3. PROTOSTELLAR EVOLUTION WITH DIFFERENT ACCRETION RATES 

There are two important timescales for evolution of accreting protostars. The first one is the accretion timescale, 

^aec^^, (1) 

over which the protostar grows by mass accretion. This is an evolutionary timescale of our calculations. The second one 
is the Kelvin-Helmholtz (KH) timescale, 

tKU ^ (2) 

over which the protostar loses energy by radiation. The balance among these timescales is crucial to the protostellar 
evolution. When > iacc, the radiative energy loss hardly affects the protostellar evolution. The evolution is controlled 
only by the effect of mass accretion. When txH < tacc, on the other hand, the evolution is mainly controlled by the 
radiative energy loss. 

Luminosity from accreting protostars has two sources: one is from the stellar interior, the other from the accretion 
shock front. In this paper, we call the former the interior luminosity L^, and the latter the accretion luminosity, 

L„ . (3) 

-ft* 

Total luminosity from the star is the sum of the interior and accretion luminosity: itot = + iacc- Note that the 
balance between txH and iacc also determines the dominant component of luminosity from accreting protostars, i.e., the 
ratio Lacc/i* is equivalent to tKn/^acc- 

During the growth of the protostar, the KH timescale significantly decreases. As a result, the ratio tKH/iacc, which is 
initially very large, falls below unity in the protostellar evolution. Below, we show that the decrease of tKH/^acc indeed 
leads to a variety of evolutionary phases of accreting protostars. 

3.1. Case with High Accretion Rate = 10-3 Mq/jv 

First, we see a protostar growing at a high accretion rate = 10-3 Mq/yi (run MD3), which is relevant to massive 
star formation. The upper panel of Figure [H shows evolution of the protostellar radius (thick solid) and the interior 
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Fig. 2. — Evolution of a protostar with the accretion rate M* = 10~^ Mg/yr (run MD3). Upper panel : The interior structure of the 
protostar. The thick soUd curve represents the protostellar radius (R*), which is the position of the accretion shock front. Convective layers 
are shown by gray-shaded area. The hatched areas indicate layers of active nuclear burning, where the energy production rate exceeds the 
steady rate Lo.st/A^* for the deuterium burning, and L,/M* for the hydrogen burning. The thin dotted curves represent the loci of mass 
coordinates ; AI = 0.1, 0.3, 1, 3, 10, and 30 Mq. Lower panel : Evolution of the mass-averaged deuterium concentration, /d,av (solid line) 
and the maximum temperature within the star Tmax (dot-dashed line). In both panels, the shaded background shows the four evolutionary 
phases ; (I) adiabatic accretion, (II) swelling, (III) Kelvin-Helmholtz contraction, and (IV) main-sequence accretion phases. 



structure as a function of the protostellar mass. The mass increases monotonically with time t, = A/*^o + M^,t, and 
can be regarded as a time coordinate. The entire evolution can be divided into the following four phases according to 
their characteristic features; (I) the adiabatic accretion (for M* < 6 Mq), (II) swelling (6 Mq < M* < 10 Mq), (III) 
KH contraction (10 Mq < < 30 Mq), and (IV) main-sequence accretion (M* > 30 Mq) phases. Below, we see the 
protostellar evolution in each phase. 



Adiabatic Accretion Phase Upper panels of Figure [3] present the evolution of entropy profile within the protostar. 
This shows that the entropy profiles at different epochs completely overlap: at a fixed mass element, the specific entropy 
remains constant. The entropy is originally generated at the accretion shock front, and embedded into the stellar interior. 
Thus, during this earliest phase, the entropy profile inside the star just traces the history of entropy at the post-shock 
point. Owing to a high value of opacity, radiative heat transport is inefficient in the interior of the protostar. This makes 
the interior luminosity small, then iKH S> tacc in this phase (see Fig. [5] bottom panel, below). Once the accreted matter 
settles into the opaque interior, the specific entropy is just conserved. 

The snapshots of the luminosity profile are also presented in lower panels of Figure [31 Near the surface is a spiky 
structure where the gradient dL/dM is very large owing to a sharp decrease in opacity there. Although a significant 
entropy loss occurs in this luminosity spike as indicated by the relation {ds/dt)M — — 1/T {dL/dM)t (from eq. (|A3P ). 
this layer is extremely thin and the time for a mass element to pass though it is so short that the entropy hardly changes. 
This can be confirmed by comparing the two timescales; the local accretion time 

Tn 

tacc,s(m) = (4) 

and local cooling time 

m Ss 

icooi,s(m) = r"^ 1 9L A~i' ' 

Jo T dm"''"' 

where m = M^ — M is the depth (in mass coordinate) of a mass shell from the accretion shock, and 6s denotes arbitrary 



6 



T. Hosokawa & K. Omukai 




2 3 
M(Mo) 




25 




150000 



aooooo 



50000 



(III) K-H contraction 



M,= 15Mo 



25 Mo 



20 Me 




1 xlO 




Fig. 3. — Radial profiles of the specifie entropy and the luminosity at different epochs for the protostar with the accretion rate M« = 
10~^ Mq/jt (run MD3) shown as functions of the mass coordinate M. For each snapshot, total stellar mass (M*) at its moment is labeled. 
The four panels correspond to the four evolutionary phases, (I) adiabatic accretion (upper left), (II) swelling (upper right), (III) Kelvin- 
Helmholtz contraction (lower left), and (IV) main-sequence accretion (lower right) phases. In the entropy profiles of the upper left panel, the 
filled circles denote the post-shock values. In the upper right panel, the open circles indicate the bottom edges of convective layers. 



small entropy change. The former, tacc.s, is the time for a thin shell of mass m at the surface to be replaced by the newly 
accreted material, while the latter, icooi,s, is the time for the same shell to lose the entropy m Ss by outward radiation. 
Comparison between these timescales in the settling layer is presented in Figure|31 where Ss — 0.1 kB/mn is adopted for 
numerical evaluation. This indicates that tacc,s is always shorter than icooi,s: the accreted material is swiftly embedded 
in the interior before losing the entropy Ss by radiation. Therefore, the adiabaticity remains valid in spite of the spike in 
luminosity profile. Note that the short tacc.s is a result of the high accretion rate. We will see that the situation changes 
for much lower accretion rates in 13.21 and 13.31 below. 

With increasing protostellar mass M^, the accretion shock strengthens and the post-shock entropy increases. Then, the 
interior entropy increases as well. This causes a gradual expansion of the stellar radius in the adiabatic accretion phase 
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Fig. 4. — Comparison between the local accretion and cooling timescales in the surface settling layer, defined by equations (|4j and 
These timescales are shown as functions of the mass coordinate measured from the accretion shock front m = M, — M. The accretion 
timescale is presented by the solid line. The cooling timescales at A/* = 1 Mq, 3 Mq, and 5 A/q are presented by the dot-dashed, dashed, 
and dotted lines, respectively. 



(Fig. upper panel). Using the typical density and pressure within a star of mass and radius i?* (e.g.. [Cox fc GiuiH 
119681) : 



and the expression for specific entropy of ideal monatomic gas 



37^, ^ P 
s = — — In 



,5/3 



const. 



(6) 



(7) 



2n VP" 

where TZ is the gas constant and /i is the mean molecular weight, the stellar radius and entropy is related as (IStahlerl 

[1983); 

i?, cx Mr^^^exp 



37^ 



(8) 



i.e., the stellar radius is larger for the higher entropy within the star. SPS86 have derived the mass-radius relation for 
protostars in the adiabatic accretion phase; 

\ 0.41 



i?, ~ 26 Rr. 



M, 



10-3 Mo/yr 



(9) 



under the condition that opacity in the radiative precursor is dominated by bound-free absorption. As this condition 
holds in our case, our calculated mass-radius relation is in a good agreement with equation Equation ^ suggests the 
large radius (> 10 Rq) of the rapidly-accreting protostar. 

The interior of the protostar remains radiative throughout this phase (Fig.[2l upper panel): all the energy transport is 
via radiation. This is in high contrast with low M■^, 10~^ Mo/yr) cases, where most of the interior becomes convective 
owing to deuterium burning for > OAMq (see Fig. [7] below). This can be attributed to a difference in the interior 
temperature. From equation ([6|), typical temperature within the star is 



II P 

r = — — 

7^ p 



G ^lAU 



(10) 



Therefore, the large stellar radius leads to the low temperature in the stellar interior. In fact, the maximum temperature 
in this phase does not exceed the threshold for deuterium burning (lower panel of Fig. [5]). 

Once the accreted matter has passed through the outermost layer with the luminosity spike, the luminosity gradient 
becomes much milder. The luminosity has a maximum at some radius: outside the luminosity maximum, the gradient 
dL/dM < 0, while dL/dM > inside. This means that heat is removed from the deep interior and absorbed in the outer 
layer. This entropy transfer, however, remains small and does not modify the entropy distribution significantly during this 
phase. Efficiency of the outward entropy transfer is related to the value of opacity. In most of the stellar interior, major 
sources of opacity are bound-bound and bound-free transitions of heavy elements, which approximately obey Kramers' 
law, K cx 

((Havashi. Hoshi fc SugimotQ.1962;.Clavton..l96&) . With the increase in stellar mass and then the interior 
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of the ZAMS stars is also presented by the thin solid line. Middle panel : Evolution of various contributions to the interior luminosity of the 
protostar. The thick dotted line represents the stellar luminosity, L» (also shown in the top panel). The total burning rate of each nuclear 
reaction is shown for the deuterium burning (Lj, coarse dotted), pp-chain (Lpp, dot-dot-dashed), and CN-cycle (Lcn, dot-dashed line). The 
horizontal coarse dotted line indicates the steady deuterium burning rate gt ■ The maximum luminosity within the star is plotted with the 
coarse dashed line. Bottom panel : Evolution of the accretion timescale tacc (dashed), and Kelvin-Helmholtz timescale (kh (dotted) for the 
same protostar. In all panels, the shaded background shows the four evolutionary phases, as in Fig. [2] 



temperature, the opacity decreases. This results in steady increase of the outward heat flux with evolution. In fact, as 
shown in the lower panels of Figure [31 the amplitude of the luminosity increases with the growth of the protostar. Also 
the maximum luminosity within the star Lmax increases as a power-law function of as seen in the middle panel of 
Figure [S] This dependence can be understood as follows: for a star with radiative stratification and Kramers' opacity, 

W 1/2 

the luminosity scales as Lrad oc M* i?* (e.g., Cox & Giuli 1968). We have confirmed that Lmax roughly obeys, 

M. \"" / R. ' 
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in our calculations. Using equations (0 and (fTTj) . we obtain Lmax oc Aff-^ for a constant accretion rate. When opacity 
becomes sufficiently small, the radiative heat transport becomes significant even in the deep interior. Consequently, the 
entropy distribution and thus stellar structure are drastically altered, which leads to the subsequent swelling phase. 



Swelling Phase At the mass of about 6 M©, the protostar begins to swell up dramatically. The stellar radius suddenly 
increases by a factor of about three and eventually exceeds 100 Rq. This swelling is caused by redistribution of entropy 
within the star. As explained above, the outward heat flow continues to increase, and ceases to be negligible. Figure 
[3] shows that the entropy distribution changes significantly with time, i.e., with M*: the entropy decreases in the deep 
interior, and increases in the outer layer. The extent of the inner entropy-losing region becomes larger and larger, and 
approaches the surface, as more of the interior becomes less opaque. 

As seen in the bottom panels of Figure [31 the peak of the luminosity, which is the boundary between the entropy- 
losing interior and the outer absorbing layer, gradually moves toward the surface with increasing height. This outward 
propagation of the luminosity peak is called the luminosity wave (SPS86). Only a thin outlying layer absorbs the entropy 
transported from the deep interior. This rapid increase of the entropy near the surface causes swelling of the radius. 
Actually, only a small fraction of mass near the surface contributes to this swelling. For example, at the maximum 
expansion of the radius i?* ~ 140 Rq at ~ 10 Mq, only 0.03% of the total mass is contained in the outer layers of 

i?* > 70 Rq. 

Note that the swelling is not due to the shell burning of deuterium. iPalla fc Stahleil (|1990( ) have presented that a 
similar swelling occuring in intermediate-mass protostars at the low accretion rate of = 10~^ Mq/ji, and concluded 
that this is caused by the shell burning of deuterium. In our case, however, deuterium burning is not important for the 
swelling. Even in the run without deuterium burning (MD3-noD), the swelling is caused by redistribution of entropy in 
the same fashion (also see § 13.31 below). During this phase, the deuterium burning occurs in the inner region (Fig. ^ 
upper panel). Despite the active deuterium burning, most of the stellar interior remains radiative. Because of the higher 
entropy and thus lower density in the star, the opacity is lower in our case whe n the active deuterium burning begins. 
This allows efficient radiative transport of the generated entropy (jStahleilll988[ ). The middle panel of Figure [5] shows 
that the maximum radiative luminosity i^ax always exceeds that by the deuterium burning Ld- This indicates a large 
capacity of radiative heat transport. Also a thin convective layer near the surface (Fig. [2]) is not brought about by the 
deuterium burning. Although this layer receives the transported entropy from inside, the higher opacity due to ionization 
prevents the outward heat flow. This causes a negative entropy gradient and then convection near the surface. 



Kelvin-Helmholtz Contraction With the approach of the luminosity wave to the surface, an increasing amount 
of energy flux escapes from the star without being absorbed inside. This results in a significant increase in the interior 
luminosity for > 5Mq (Fig. (5]) and thus the shorter KH timescale ^kh {oc 1/L*) of the star. The arrival of 
the luminosity wave at the stellar surface marks the moment that all parts of the star begin to lose heat. Almost 
simultaneously, ^kh becomes as short as the accretion timescale tacc (Fig-E])- While maintaining the virial equilibrium 
against the radiative energy loss, the protostar turns to contraction. This is the so-called KH contraction phase. The 
epoch of the radial turn-around is roughly estimated by equating tacc and tKH.imax = GAf^/i?*Linax, where imax is given 
by equation pT|) . 

/ . V 2/9 

M™ = 10.8 Mo (^^^^,3^ j . (12) 

In the above, we have omitted a weak dependence on i?^ and just substituted i?* ^ 10 Rq for simplicity. During 
the contraction of the protostar, the KH timescale remains slightly shorter than the accretion timescale. Since the ratio 
ixH/^accj or equivalently Ljf/L„iax) is iess than unity, the interior luminosity exceeds the accretion luminosity and becomes 
the dominant source of the total luminosity of the protostar (upper panel of Fig. [5]). 

Since the stellar interior remains radiative in spite of the deuterium burning, the accreted deuterium is not transported 
to the deep interior by convection. Deuterium in the deep interior is thus soon consumed. Thereafter the deuterium 
burning continues in a shell- like outer layer (Fig. [21 upper panel). This surface burning proceeds at a rate slightly 
exceeding the steady-burning one fe.g.. IPalla fc Stahleri[l99l . 

. . MJ, = 1500 Lq ( f'* , \ ( P/^] \ . (13) 
° \^ 10-3 Mg/yry V 2-5 X 10-5 y ' ^ 

where Su is the energy available from the deuterium burning per unit gas mass. Consequently, the mass-averaged deuterium 
concentration /d.av significantly decreases during the KH contraction (Fig. [2l lower panel). 



Arrival at a Main-sequence Star Phase With the KH contraction, the interior temperature increases gradually. 
At Af, ~ 20 Mq the maximum temperature reaches lO^K and the nuclear fusion of hydrogen begins (Fig. [5]). Although 
the hydrogen burning is initially dominated by the pp-chain reactions, energy generation by the CN-cycle reactions 
immediately overcome this. The energy production by the CN-cycle compensates radiative loss from the surface at 
A'f, ~ 30 Mq (Fig. O middle), where the KH contraction terminates. A convective core emerges owing to the rapid 



10 



T. Hosokawa & K. Omukai 



1000 




o 



100 



10 
10 



0.1 



U = 10-3 M®/yr 








: accretion shock 






: photosphere 
















■ ■ ■ ■ 1 


/ 








/ . 
/ / 








/ / 








/ / 




Tprec / 10 




1 / 
1/ 
1/ 

/ Tsh,4 




H . 






\ 


^ 1 


Tph,4 




I 
1 
I 


1 
1 






1 

• ■ 


1 
1 





10 

logM,(Me) 



Fig. 6. — Upper panel : Positions of the accretion shock front (sohd hne), and photosphere (dotted hne) of a growing protostar with 
M» = 10~^ Mq/jt (run MD3). The layer between the accretion shock front and photosphere corresponds to the radiative precursor. Lower 
panel : the optical depth within the radiative precursor Tree, and effective temperatures at the photosphere Tp^ and accretion shock front 
= {Lf /iirR^a-)^^^ . The temperatures are normalized as T4 = (T/lff^K). The shaded background shows the four evolutionary phases, as 
in Fig. [2] 



entropy generation near the center. The stellar radius increase after that, obeying the mass-radius relation of main- 
sequence stars. 

Evolution of the Radiative Precursor Figure [S] shows the evolution of protostellar (i?,) and photospheric (-Rph) 
radii. Except in a short duration, the photosphere is formed outside the stellar surface: the radiative precursor persists 
throughout most of the evolution. In the adiabatic accretion phase, the photospheric radius is slightly outside the 
protostellar radius and increases gradually with the relation i?ph — 1.4i?*, which is derived analytically by SPS86. 
Although the precursor temporarily disappears in the swelling phase around ~ IOMq, it emerges again in the 
subsequent KH contraction phase. In the KH contraction phase, the precursor extends spatially and becomes more 
optically thick. Due to the extreme temperature sensitivity of H~ bound-free absorption opacity, the region in the 
accreting envelope with temperature > 6000K becomes optically thick. Thus, the photospheric temperature Tph remains 
at 6000 — 7000K throughout the evolution (Fig. [6l lower panel). Since the photospheric radius is related to the total 
luminosity by 

i?ph = (itot/4^fTrp\)'/', (14) 

for the constant Tph, the rapid increase in the luminosity in the contraction phase causes the expansion of the photosphere. 

3.2. Case with Low Accretion Rate = 10^^ Mq/jt 

Next, we revisit the evolution of an accreting protostar under the lower accretion rate of A/* = 10^^ M^g/yr (run MD5). 
Although this case has been extensively studied by previous authors, a brief presentation should be helpful to underline 
the effects of different accretion rates on protostellar evolution. Different properties of the protostar even at the same 
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Fig. 7. — Same as Fig. [2] but for the lower accretion rate of M» = 10~^ Mq/jt (run MD5). In the lower panel, deuterium concentration 
in the convective layer cv is also presented. In both upper and lower panels, the shaded background shows the four evolutionary phases ; 
(I) convection, (II) swelling, (III) Kelvin-Helmholtz contraction, and (IV) main sequence accretion phases 



accretion rate owing to updates from the previous works, e.g, initial models and opacity tables, will also be described. 
For a thorough comparison between our and some previous calculations, see appendix iBl 

The lower accretion rate means the longer accretion timescale tacc cx; even at the same protostellar mass M^^ 

the evolutionary timescale is longer in the case of low M* . Therefore, the protostar has ample time to lose heat before 
gaining more mass. This results in the lower entropy and thus a smaller radius at the same protostellar mass in the low 
M* case as shown in the upper panel of Figure [71 Although with the smaller value, the overall evolutionary features of 
the high and low radii are similar. 

Again, the protostellar evolution can be divided into four characteristic stages, i.e., (I) convection (M* < 3 Mq), 
(II) swelling (3 Mg < M* < 4 Mq), (III) KH contraction (4 M© < M, < 7 Mq), and (IV) main-sequence accretion 
{M^, > 7 Mq) phases. Note that the first phase is now the convection phase instead of the adiabatic accretion. In the low 
M* case, the deuterium burning drastically affects the protostellar evolution in early phases fe.g.. IStahledllOSSl ). Since 
the evolution in the phases (III) and (IV) is similar to the previous high-M* case 13. ip . we emphasize the two early 
phases (I) and (II) in the following. 

Convection Phase In the low M* case, the D burning begins in an earlier phase and has a more profound effect on 
the evolution than in the high-M* case (run MD3). As a result of the smaller radius, the temperature, as well as density, 
is higher in the low M, case (cf. lower panels of Fig. [2] and Fig. [7|). The maximum temperature reaches 10^ K as early 
as M* ~ 0.2 Mq and active deuterium burning begins subsequently (Fig. [7]), while in the high-M* case the D ignition 
does not occur until M* ~ 6 Mq (Fig. [J). Unlike in the high-M* case, the deuterium burning makes most of the interior 
convective. This is a result of high density and thus high opacity at the ignition of deuterium burning, occurring at a 
fixed temperature of ~ 10^ K. Since high opacity hinders the radiative heat transport, entropy generated by deuterium 
burning cannot be transported radiatively. This necessarily gives rise to convection. 

Figure [7] shows that the extent of the convective layers experiences a complicated history. At M* ~ 0.3 Mq, a convective 
layer appears for the first time slightly off-center. Inside the convective layer, entropy is homogenized and its value increases 
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Fig. 8. — Same as Fig.[8]but for the lower accretion rate of Af« = 10~^ Mq/jx (run MD5). 
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owing to the D burning (Fig. [51 1). With this increased entropy, the outer layer is graduaUy incorporated into the convective 
region. Fresh deuterium is supphed from the newly incorporated layer and spreads over the convective region, and finally 
is consumed by nuclear burning at the bottom of the region. However, the expansion of the convective region cannot be 
maintained if the supply of deuterium becomes insufficient. Without the fuel, the deuterium concentration in the convective 
region /d,cv drops and the nuclear burning ceases: the convective region disappears. This occurs at Af* ~ 0.6 Mq (Fig. [71 
upper). At ~ 0.7 Mq, another episode of deuterium burning begins just outside the first convective region, where 
fresh deuterium still remains. The second convective region develops there and immediately reaches the stellar surface. 
Although this complicated evolution of the convective regions is different from that reported in SSTSOb, this difference 
can be attributable to high sensitivity of the evolution on adopted initial deuterium abundance. We discuss this in greater 
details in Appendix lB.il 

Since the energy generation by D burning is very sensitive to the temperature, temperature remains constant at 
~ 1.5 X lO^K during the active deuterium burning (Fig. [71 lower panel). This is the so-called thermostat effect (e.g., 
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Fig. 9. — Same as Fig.[8]but for switching off deuterium burning in the calculation (run MD5-noD). Only profiles in the early adiabatic 
accretion phase and swelling phase are presented. 



Stabler 1988). If this works, M^jR^ oc T ~ const, from equation (fTO|) : the protostellar radius i?, increases as oc M*. This 
linear increase in i?* with can be observed in the range O.TMq < < IMq (Fig. [71 upper panel). 

When the deuterium concentration in the convective region drops to a few percent, the thermostat effect ceases to work. 
The maximum temperature starts increasing again and the radial expansion slows down after that. Owing to the high 
temperature, opacity decreases in the deep interior with increasing stellar mass. This allows efficient radiative energy 
transport and renders a large portion of the interior radiative. The radiative core gradually extends to the outer layer, and 
the convective layer moves outward. Since fresh deuterium is no longer supplied to the radiative interior by the mixing, 
deuterium is soon exhausted there. The expansion of the radiative interior occurs in a later phase in the calculation of 
PS91. This difference can be attributed to different initial models: whereas PS91 assumed a fully convective star of 1 Mq 
as an initial model, we started calculation from the initial mass M^^q = 0.01 Mq. At M* = 1 Mq, the protostar is not 
fully convective in our case. The outer layers of M > 0.25 Mq are convective, where the D burning occurs at the bottom, 
while the remaining inner core is radiative at this moment. We confirmed that, starting from the same initial model, the 
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Fig. 10. — Same as Fig. |4] but for the case with the lower accretion rate of M, = 10~^ Mq/yv and without deuterium burning (run 
MD5-noD). The local cooling timescales given by equation (|5j at Af, = 1.0 (dot-dashed), 1.5 (dashed), and 2.0 Mq (dotted line) are 
presented. 
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Fig. 11. — Same as the upper panel of Fig.[7]but for the case without the deuterium burning ([D/H] = 0, run MD5-noD). The dot-dashed 
hne represents the mass-radius relation for the case with the fiducial deuterium abundance, [D/H] = 2.5 X 10~^ (run MD5, also see Fig. [7]l. 



evolution becomes similar to that of PS91 as discussed in Appendix IB. 21 

Because of the early ignition of deuterium burning, entropy and luminosity profiles in this phase look rather different 
from those in the high case (see Fig. [3] and |5]). For easier comparison, let us see the case without deuterium burning, 
where the similar profiles are in fact reproduced (Fig.[21). However, one important difference exists, the spiky structure in 
the entropy profiles near the surface. This shows significant entropy loss near the surface in the low M* case, which leads 
to the low entropy within the star. The reason can be seen clearly in Figure [TOl where the comparison between the local 
accretion timescale tacc,s and cooling timescale ^cooi.s (eq.s|3]and[ni) is presented. In contrast to the the case of high M*, 
in most of the settling layer, tacc,s is shorter than icooi,s: accreted matter loses entropy radiatively near the surface before 
settling into the adiabatic interior. 



Swelling Phase In the short period of 3 M© < Af, < 4 Af©, the stellar radius jumps up by a factor of two. At the 
same time, the whole interior becomes radiative (Fig. [71 upper). Although deuterium burning continues near the surface, 
it does not drive convection any more. 

The cause of this swelling has been attributed to the shell burning of deuterium, which also occurs in our calculation 
(Palla & Stabler 1990). However, we conclude that the main driver of the swelling is not the D shell burning, but the 
outward entropy transportation as in the case of the high accretion rate. The reasons are as follows. First, evolutionary 
features of entropy and luminosity profiles in this phase are similar to those in the swelling phase in the high case 
(run MD3; cf. FigslSllH]): in both cases, entropy decreases in the deep interior and increases significantly near the stellar 
surface; simultaneously, the peak of the luminosity distribution moves to the surface as the luminosity wave. Second, even 
in the run without deuterium burning (run MD5-noD), this swelling occurs as shown in Figure [TT] the protostar swells 
up at A//* ~ 3 Mq. 

The interior luminosity jumps up at Af, ~ 0.7 Mq with the arrival of the convective region at the surface as a result 
of high efficiency of convective heat transfer (Fig. 1121 upper panel). Although this rise in decreases tKH, ^kh remains 
much longer than t^cc (Fig. 1121 lower panel) . It is not until another sudden rise in at ~ 3 — AMq due to the arrival of 
the luminosity wave that tKH becomes as short as iacc- The maximum radius is attained at this epoch since the swelling 
is also caused by the luminosity wave, as seen above. Then, the analytical argument leading to equation (|12p also holds 
in this case: for Af, = 10^^ Mq/ji, this leads to the maximum radius at Af,^rmax — 3.9 Mq, which agrees with our 
numerical result. 

Simultaneous returning to the radiative structure with the swelling has been explained bv lPalla &: Stahleil ()1990D . The 
luminosity from the steady deuterium burning is comparable to the accretion luminosity by chance: 

Lo.st = M^Sr, ~ M^— — = Lace- (15) 

The epoch of the radial turn-around obeys equation (|12p . which is derived by equating tacc and iKH,imax — GAf^/i?*Lmax- 
This means that Lmax becomes comparable to Lace and thus to Lo.st in the swelling phase. Hence, the radiative luminosity 
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Fig. 12. — Same as Fig. [5] but for the lower accretion rate of Af» = 10~^ MQ/yr (run MD5). In eacli panel, the stellar luminosity L, 
(Top ■panel, thin dotted line), maximum luminosity within the star Lmax {Middle panel, thin coarse dashed line), and Kelvin-Helmholtz time 
tKH {Bottom panel, thin dotted line) for the case without deuterium burning (run MD5-noD) are presented for comparison. In all panels, the 
shaded background shows the four characteristic phases, as in Fig. [71 



now has a capacity to transport the energy generated by deuterium burning: Lmax ^ iD,st- The deuterium burning no 
longer gives rise to any convective layer after the swelling. 

Later Phases Subsequently, the protostar enters the KH contraction phase. The evolution thereafter is similar to the 
previous high-M* case. The maximum temperature within the star increases with the contraction, and hydrogen burning 
begins at M* ~ 7 M©. 

Evolution of the Radiative Precursor Figure [T51 shows that the accretion flow is optically thin except in the early 
phase of M* < 1 AIq at the low accretion rate of 10~^ Mq/yi. With the expansion at the active deuterium burning, the 
stellar surface emerges above the photosphere at ~ 1 Mq and the radiative precursor disappears thereafter. Without 
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Fig. 13. — Same as Fig. |6] but for the lower accretion rate of Af* = 10~^ Mq/jt (run MD5). The shaded background shows the four 
characteristic phases, as in Fig. [7] 

the precursor, the steUar surface is directly visible, with effective temperature rising remarkably in the KH contraction 
phase, i.e., for M# > AMq. This is in stark contrast to the case of high accretion rate 10~^ M©/yr, where the effective 
temperature remains several thousand K throughout the evolution (Fig. [S]) . 

3.3. Dependence on Mass Accretion Rates 

Here, we summarize the effects of the different accretion rates on the protostellar evolution. The mass-radius relations 
of protostars with accretion rates M* = 10~^, 10~^, 10~^, and lO~^M0/yr are shown in the upper panel of Figure [Ml 

First, the protostellar radius is larger for the higher accretion rate at the same protostellar mass, as discussed in ij 13.11 
For example, at = IM©, the radius i?, ~ 25i?0 for 10~'^Mo/yr, while it is only 2.5Rq for lO~^Af0/yr. For easier 
comparison, the mass-radius relations in the no-deuterium cases are also shown by thin lines in Figure [HI The deuterium 
burning enhances the stellar radius, in particular, in low cases (see in this section below). Nevertheless, the trend of 
larger i?, for higher A/* remains valid even without the deuterium burning. For simplicity, we do not include the effect of 
deuterium burning in the following argument. The origin of this trend is higher specific entropy for higher accretion rate. 
Recall that the stellar radius is larger with the high entropy in the stellar interior as shown by equation ([8]) . The entropy 
distributions at = IMq for the no-D runs are shown in Figure jTHl which indeed demonstrates that the entropy is 
higher for the higher cases. The interior entropy is set in the post-shock settling layer where radiative cooling can 
reduce the entropy from the initial post-shock value as observed in spiky structures near the surface (e.g., the case of 
M* < 10~^ Mq/ji in Fig. [T5)) . For efficient radiative cooling, the local cooling time icooi,s in the settling layer must be 
shorter than the local accretion time tacc,s, as explained in i) 13.11 and 13.21 Figure [TBI shows that the ratio icooi,s/iacc,s is 
smaller for lower and becomes less than unity for < 10"^ Mq/yi. Thus the accreted gas cools radiatively in the 
settling layer in low (< lO~^M0/yr) cases, while in the higher M, cases, the gas is embedded into the stellar interior 
without losing the post-shock entropy. The difference in entropy among high Af, (> lO^^Afo/yr) cases comes from the 
higher post-shock entropy for higher Af*. In those cases, the protostar is enshrouded with the radiative precursor, whose 
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Fig. 14. — Evolution of the protostellar radii (upper panel) and the maximum temperatures within the stars (lower panel) for different 
mass accretion rates. The cases of Af* = 10"*^ (run MD6; solid), lO'^ (MD5; dashed), 10"* (MD4; dot-dashed), and IQ-^Me/yr (MD3; 
dotted line) are depicted. Also shown by the thin lines are the runs without deuterium burning ("noD" runs) for the same accretion rates. 
In both panels, all the curves finally converge to single lines, which is the relations for the main-sequence stars. 



optical depth is larger for the higher accretion rate (c.f..FigsI^ and \W[ bottom). Consequently, more heat is trapped in 
the precursor, which leads to the higher post-shock entropy. 

Figure [Til shows not only that not only the swelling of the stellar radius occurs even in the "no-D" runs, but also that 
their timings are almost the same as in the cases with D. This supports our view that the cause of the swelling is outward 
entropy transport by the luminosity wave rather than the deuterium shell-burning I3.1l and l3.2p . Also, in all cases, the 
epoch of the swelling obeys equation (|12p . which is derived by equating the accretion time tacc and the KH time iKH,imax- 

Another clear tendency is that the deuterium burning begins at higher stellar mass for higher accretion rate. Similarly, 
the onset of hydrogen burning and subsequent arrival at the ZAMS are postponed until higher protostellar mass for the 
higher accretion rate. For example, the ZAMS arrival is at M, ~ AMq {AOMq) for 10^^ (10^'^, respectively) Mq/jy. 
These delayed nuclear ignitions reflect the lower temperature within the protostar for the higher at the same M* 
(Fig. 1141 lower panel). Since i?, is larger with higher at a given Af,, typical temperature is lower within the star. 
Substituting numerical values in equation (jlOp . we obtain 

which is a good approximation for our numerical results. 

Finally, we consider the significance of deuterium burning for protostellar evolution. Comparing with "no-D" runs, we 
find that the effect of deuterium burning not only appears later, but also becomes weaker with increasing accretion rates. 
At the low accretion rates M* < 10~^ MQ/yr, convection spreads into a wide portion of the star owing to the vigorous 
deuterium burning. In addition, the thermostat effect works and the stellar radius increases linearly during this period. 
At the higher accretion rates Af* > 10"'' Mq/jt, on the other hand, deuterium burning affects the protostellar evolution 
only slightly. Convective regions do not appear even after the beginning of deuterium burning for A/* = lO~^Af0/yr 
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Fig. 15. — Distributions of the specific entropy for the case without deuterium ("no-D" runs in Table 1) at the epoch of Af» = 1 Mq for 
different accretion rates. Histories of the post-shock values arc also plotted by the dotted lines for the cases with Ah = 10~^ and 10^^ Mg/yr. 
For the cases with M* = 10""* and 10"'' Mg/yr, those curves completely overlap with the entropy distributions at Mt = 1 Mq. 
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Fig. 16. — Evolution of the minimum value of the ratio tcool,s("^)/tacc,s (m) in the surface settling layer. The no-D cases are shown for 
Alt = 10~® (solid), 10~^ (dashed), 10~* (dot-dashed), and 10~^ Mg/yr (dotted). For the ratio icool s(™)/tacc,s (™) exceeding unity, the 
accreted material settles into the interior without losing its post-shock entropy radiatively. 



(§1311). 

A key quantity for understanding this variation is the maximum radiative luminosity at the ignition of deuterium 
burning. The maximum radiative luminosity within the star Lmax at deuterium ignition can be evaluated by eliminating 
Rt and in equation ^ with (H]), HI]), and setting T^ax ^ 10^ K: 

^--^»^^^( iO-3M,/yr ) • 

On the other hand, the luminosity from deuterium burning is (eq. [T3|) Ln.st ~ 10'^ Lq {M^,/10~^ Mo/yr). For accretion 
rate > IQ-^ Mg/yr, L 

max > ^max.D- radiation is able to carry away all the entropy generated by the deuterium 
burning. Therefore, the protostellar interior remains radiative in the high M* cases (e.g., run MD3, see § 13. ip . while 
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Fig. 17. — The protostellar mass-radius relations for different metallicities. The cases of the solar (dotted; runs MD5 and MD3) and zero 
(solid; MD5-zO and MD3-zO) metallicities are presented for the accretion rates M* = 10~^ and lO~^M0/yr. 



convection is excited for the lower (e.g., run MD5, see 13. 2p . 

3.4. Comparison with Primordial Star Formation 

Besides the massive star formation in the local universe, high accretion rates have been invoked in the primordial star 
formation in the ear ly universe. Without important metal coolants, the temperature remains at several h undred K in the 
primordial gas fe.g.. lSaslaw fc Zipovl[l967HMatsuda. Sato fc Takedalll969HPalla. Salpeter k Staliiei1ll983D . The accretion 
rate, which is given by 

^3 / y \3/2 

^ = lO-^Afg/yr — - , (18) 



G V500K, 

where Cs and T are the sound speed and temperature, respectively, in star-forming clouds, becomes as high as 10^^ Mq/jt 
even without turbulence (SPS86). As a result of similar accretion rates, protostellar evolution of the primordial stars well 
resembles with those of the high M* cases studied in this paper (SPS86 ; Omukai & Palla 2001,2003). The mass-radius 
relations for the primordial {Z ~ 0) protostars are presented in Figure [T71 along with their solar-metallicity {Z = 0.02) 
counterparts for comparison. This indicates that basic evolutionary features are similar despite differences in metallicity: 
as in the present-day protostars, primordial protostars pass through the same four characteristic phases (see § 13.11 and 
13. 2p . The following two differences, however, still exist. First, the swelling of the protostar and subsequent transition to 
the KH contraction occur earlier in the Z = cases. For example, for M^, = lO~^M0/yr, the maximum radius is attained 
and the KH contraction begins at IOMq in the Z = 0.02 case, while it is 7Mq in the Z = case. Recall that the transition 
to those phases is caused by the decrease in opacity with increasing Af, (§ 13. ip . Owing to the lack of heavy elements, 
the opacity is lower in the primordial gas at the same density and temperature. Efhcient radiative heat transport begins 
earlier, which results in the earlier transitions in the evolutionary phases. Second, the stellar radius in the main-sequence 
phase is smaller as a result of higher interior temperature (~ lO^K) for the primordial stars. Due to the lack of C and 
N, sufficient energy production by the CN cycle is not available at a temperature similar to that of the solar metallicity 
case (~ lO^K). The KH contraction is not halted until the tempera ture reaches c± lO^K, where a small amount of carbon 
produced by the He burning enables the operation of the CN cycle (|Ezer fc Cameronl[l97TI : lOmukai fc Pallall2003h . bmce 
more contraction is needed, the arrival at the ZAMS can be later in the Z = case, despite the earlier beginning of the 
KH contraction: for example, it is 30Mo {50Mq) in the Zq (Z = 0, respectively) case. 

3.5. Radiation Pressure Barrier on the Accreting Gas Envelope 

The higher the mass accretion rate, the more massive the protostar grows before the arrival at the ZAMS phase (§ 13.11 
and l3.3p . However, this trend does not continue unlimitedly: for a rate exceeding a threshold value, stars have an upper 
mass limit set by the radiation pressure onto the radiative precursor. 

The evolution of the protostellar radius is shown in Figure [TH] (upper panel) for cases with very high accretion rates 
> 10~^ Mq/jt. Up to a certain moment in the KH contraction phase, the evolution in all cases proceeds in a qualitatively 
similar way although the radius is larger and the onset of swelling and the KH contraction are later for higher Af, 
cases. However, for accretion rates M* > 3 x 10~^MQ/yr, the contraction is halted at ~ 40Mq, followed by abrupt 
expansion at ~ 50Mq. In such cases, further steady accretion onto the star is not possible owing to strong radiative 
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Fig. 18. — Upper panel : Evolution of the protostellar radii for cases at the very high accretion rates of 10 ^ (soUd, run MD3), 3 X 10 ^ 
(dashed, MD3x3), 4 X IQ-^ (dot-dashed, MD4x3), and 6 X IQ-^ Mg/yr (dotted, MD6x3). For comparison, we also plot the case of a 
primordial {Z = 0) protostar with 3 X 10~^ Mq/jt (thin dashed, MD3x3-zO). Lower panel : Evolution of a ratio between the total luminosity 
of protostars and Eddington luminosity for the same cases. 



pressure exerted onto t he ra diative precursor. This phenomenon has been found for the primordial star formation by 
lOmukai fc Palial (|200lL l2003f) . They found that abrupt expansion occurs when the total luminosity from a protostar, 
Ltot = i* + iacc = + GM^M^,/ R^, becomes close to the Eddington luminosity ^Edd- The strong radiation pressure 
decelerates the accretion flow and then reduces the ram pressure onto the stellar surface. With reduced inward ram 
pressure, an outward thrust by the interior radiation blows a thin surface layer away. 

The critical accretion rate Af* is defined as the rate above which the abrupt stellar expansion occu rs in the KH 
contr action phase. By using the condition that Ltot reaches i^Edd before the protostar reaches the ZAMS, (|Omukai fc Palial 
l2003f ) derived the critical rate as follows: 

. GMzAMsM^.cr . 

-f^ZAMS H J5 = -t^Edd, (iyj 

^ZAMS 



that is, 

-7 _ 47rcjT;zAMS L 

-'W*,cr — I 1 ; 

Kosc V ^Edd 

where k^sc is the electron-scattering opacity, and the suffix "ZAMS" means quantities of the ZAMS stars. Although 
the critical rate M^^ci by equation ((20|) depends on the stellar mass, the dependence is very weak and M^ cr — 4 x 
10~'^AfQ/yr for the primordial stars. Since the ZAMS radius increases with metallicity (see 13. 4p . M*.cr is expected to 
increase with metallicity. For example, for Z — 10~^, M*^cr = 9 x lO~'^Af0/yr by equation ([20l) . and the evolutionary 
calculation demonstrated that abrupt stellar expansion occurs in fact for higher accretion rates (|Omukai fc Palla,20()3l ). 
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Our calculations, however, indicate that the critical accretion rate in the solar metallicity case is ~ 3 x lO^'^M0/yr, 
slightly smaller than that of primordial protostars. This apparent contradiction comes from the fact that, in the solar 
metallicity case, opacity in the accreting envelope is higher than that of the electron scattering one Kos postulated in the 
above derivation for the critical accretion rate. In fact, the violent stellar expansion is found to take place with a total 
luminosity of about a half of the Eddington luminosity, as shown in the lower panel of Figure [THl A marginal case is 
Run MD3x3, where the contraction is about to be reversed, but the star barely reaches the ZAMS after some loitering. 
In this case, the Eddington ratio /Edd = itot/^Edd slightly exceeds 0.5 at ~ 45 Mq. For comparison, the case of 
the primordial star is shown for the same accretion rate, = 3 x 10~^ Mq/ji, which is slightly lower than the critical 
accretion rate (run MD3x3-zO). Because of the earlier transition to the KH contraction phase, the Eddington ratio also 
begins to increase earlier than the present-day counterpart (run MD3x3). Despite high luminosity almost reaching the 
Eddington limit, the KH contraction is not reversed and the protostar can grow further. Taking into account the higher 
opacity in the envelope, we add the Eddington factor of 0.5 to equation 

GA/zAMsM*.cr ^, ^ ^91^ 

^ZAMS H 5 — U.Oi^Edd (.^tj 

-KZAMS 

for present-day protostars. The critical accretion rate becomes M*,cr — 3 x lO~^M0/yr, which agrees well with our 
numerical results. Although this critical rate is a function of the stellar mass, its dependence is weak as in the primordial 
case. 

Our results suggest an upper mass limit of PMS stars at ~ 60 Mq , which corresponds to the mass where the protostar 
reaches the ZAMS at the critical accretion rate of Ma- 3 x 10~^ Mq/yi. More massive PMS stars cannot be formed by 
steady mass accretion. With higher accretion rates, the steady accretion is halted by the abrupt expansion at ~ 50 Mq. 
What happens after the abrupt expansion is speculated upon in the next section. 

4. UPPER LIMIT ON THE STELLAR MASS BY RADIATION FEEDBACK 

As long as the central star is not so massive and its feedback on the accreting envelope is weak, the accretion proceeds at 
a rate set by the condition of the parent star-forming core. Our assumption concerning the accretion at a pre-determined 
rate remains valid. However, when the central star grows massive enough and its feedback, e.g. radiative pressure on 
the dusty outer envelope, thermal pressure in the HII region, etc., can no longer be neglected, the accretion rate is now 
regulated by the stellar feedback. The mass accretion is eventually terminated and the final stellar mass is set at this 
moment. 

Among possible types of feedback, the radiation pressure on the outer dusty envelope is considered to be the most 
important in terminating the protostellar accretion (e.g., WC87). Since no such outer envelope has been included in our 
calculation, we evaluate here the epoch of accretion termination by stellar feedback with the assumption that the accretion 
rate is roughly constant up to this moment. 

Most of the stellar radiation, emitted in the optical or UV range, is absorbed in a thin innermost layer of the dusty 
envelope, i.e., the dust destruction front, and the outward momentum of photons is transfered to the flow (see Fig. [T]). 
Overcoming this radiation pressure barrier is a necessary condition for continuing accretion. The condition that the 
outward radiation pressure falls below the inward ram pressure of the accretion flow is expressed as 

where i?d is the radius of the dust destruction front, and p and u are the flow density and velocity there. The limit on 
the stellar mass by this condition is shown in Figure [T9l as a function of the mass accretion rate (i.e, "Pram limit"). Here, 
we have assumed that the flow is in the free-fall. In Figure [TOl the epochs of the maximum swelling and the arrival at the 
ZAMS are also presented from our results. Except in cases with very high accretion rates (> 3 x 10""^ Mq/yt; it 13. 5p . the 
radiation pressure limit is attained after the stars reach the ZAMS. This verifies previous authors' treatment in assuming 
that the central star is already in the ZAMS phase when the feedback becomes important (e.g., WC87). 

In reality, the accretion flow is decelerated before reaching the dust destruction front. The radiation absorbed at the 
dust destruction front is re-emitted in infrared wavelengths and this re-emitted radiation imparts pressure onto the outer 
flow. This radiation pressure exceeds the gravitational pull of the star if the luminosity exceeds (McKee & Ostriker 2007) : 

incGM^ f M^\ 

iEdd,d = ^ 16OOL0 — ^ , (23) 

Kd.IR \Mq J 

which is the "Eddington luminosity" for dusty gas, calculated by using the dust opacity for the far-infrared light Kd.iR — 
8 cm^/g instead of the electron-scattering opaci ty. The adopted valu e 8 cm^/g is the Rosseland mean opacity at T ~ 600 K, 
where the opacity takes the maximum value (jPoUack et al.l [199^ . In the outermost region of the dust cocoon, where 
T < 600 K and Kd,iR < 8 cm^/g, the deceleration effect is somewhat milder. This effect becomes significant when the flow 
reaches the inner warm region. As shown by the light-shaded region in Figure [HI the deceleration of dusty flow occurs 
(Ltot > iEdd,d) in a wide range of and M*. In this case, the upper mass limit by the radiation pressure ("Pram limit" 
in Fig. [in]), which is derived by assuming the free-fall flow, becomes lower: the boundary of the "Pram limit" region shifts 
to the left to some extent. To evaluate this effect more quantitatively and then to identify how mass accretion is finally 
terminated, further studies will be needed. It has been shown that the deceleration effect of infrared light is diminished by 



22 



T. Hosokawa & K. Omukai 




Fig. 19. — Limits on the stellar mass for each accretion rate by various feedback effects. The three horizontal arrows represent protostellar 
evolution at the accretion rate of (i) M* = IQ-^ (run MD5), (ii) 10"^ (run MD3), and (iii) 6 X 10"^ Mq/yx (run MD6x3). Also shown 
are the stellar masses where the protostellar radii reach the maximum (dashed line) and where the protostars arrive at the ZAMS (dotted 
line). The region shaded by light gray represents where the total luminosity of the protostar, Ltot = L* + Lace exceeds the dust Eddington 
luminosity LEdd.di given by equation II23I I. The lower-right hatched region represents where the ram pressure of the free-fall flow is insufficient 
to overcome the radiation pressure barrier at the dust destruction front. The upper-right hatched region denotes the similar prohibited region 
but for the radiation pressure acting on a gas envelopes just around the protostar. The extra arrow (iv) indicates a possible path along which 
the protostar will evolve after the gas envelope is disrupted by the radiation pressure (see text). 



reduced dust opacity (WC87), non-spherical accretion (e.g.. lNakanolll989t lYorke fc SonnhalteHl2003 ). and other various 
processes (e.g., see McKee & Ostriker 2007 for a recent review). Note that even with the reduced dust opacity, the "Pram 
hmit" still remains as long as the dusty envelope is thick enough for optical/UV (not infrared) light from the star. 

As seen in § 13.51 for a very high accretion rate of M* > 3 x 10"'^ Mq/jt, the radiation pressure acting on the gas 
envelope terminates the steady mass accretion. Figure [19] shows that this limit ( "quasi-I/Edd limit" ) is more severe than 
the radiation pressure limit for > 3 x 10~^ Mq/ji although the deceleration of the dusty flow is already important. 
Note that WC87's result on this quasi- Eddington limit differs from ours (see their Fig. 5). This discrepancy comes from 
their assumption that the forming star is already in the ZAMS, which is not valid in that case. On the other hand, 
we have consistently calculated the evolution until this limit is reached and the star begins abrupt expansion. For such 
a high accretion rate, what happens when the steady accretion becomes impossible is not clear in realistic situations. 
The accreting envelope may be completely blown away and growth of the protostar may be terminated at this point. 
Without accretion, the protostar quickly relaxes to a ZAMS star. Another possibility is that since the strong radiation 
force originates in the accretion luminosity itself, the accretion may continue in a self-regulated way by reducing its rate, 
for example, by an outflow. This process may proceed in a sporadic way: the protostar continues to grow alternately 
repeating eruption and accretion. If this is the case, the protostar evolves along the trajectory (iv) in Figure 1191 which 
runs along the boundary of the prohibited region and the maximum stellar mass is about 250 Mq set by the radiation 
pressure feedback at the critical accretion rate M*_cr — 3 x 10"'^ Mq/jx (Fig. [T9)) . This is a somewh at higher value than 
the likely upper mass cut-off of the Galactic initial mass function ~ 150 Mq (e.g.. IWeidner fc Kroupa. 2004, : .Figcr. 20051) ■ 



Evolution of Massive Protostars 



23 




4.5 4 3.5 3 2.5 

logT.,eff(K) 

Fig. 20. — The stellar birthlines for different accretion rates. The cases for accretion rates 10~® Mirj/yr (solid; run MD6), 10~^ (dashed; 
MD5), 10"* (dot-dashed; MD4), and IQ-^ Mq/jv (dotted; MD3) are presented. Each track shows the evolution from the initial model, and 
filled circles on the track represent points of the stellar mass M<, = 1, 3, 5, 9 , and 20Mp, from the lower right in this order. The thick broken 
line represents the positions for zero-age main-sequence from lSchaller et al.l II1992I) . The open squares along the line denote the positions for 
the same masses as above. The dotted straight lines indicate the loci for the constant stellar radius of IRq , lOiiQ , and lOOi?0 . 



This disagreement might be due to the deceleration of the outer dusty envelope, which is not included in this analysis. 

5. OBSERVATIONAL SIGNATURES OF HIGH ACCRETION RATES 

Massive protostars with high accretion rates have some characteristic observational features. In this section, we explore 
possibilities of identifying them observationally and of verifying the postulated high accretion rates. Although protostars 
still in the main accretion phase are hard to observe directly, pre-main sequence (PMS) stars just after this phase are 
optically visible. The distribution of such PMS stars in the Hertzprung-Russell (HR) diagram retains a footprint of the 
previous accretion phase. The birthline is the initial loci of the PMS stars on the HR diagram as a function of the 
stellar mass. The PMS stars move down-left toward the ZAMS line on the diagram thereafter. It has been shown that 
the birthline for M* ~ 10~^ Mq/jy successfully traces the upper envelope of observed distrib utions of low-mass and 
intermediate-mass PMS stars in the HR diagram (iPalla fc Stahled[T99^ : lNorberg fc Maedeill2000D . In Figure HOI we show 
the birthlines for various accretion rates. These tracks are calculated from the interior luminosity and the stellar radius 
R^, at the end of the accretion, i.e., when the protostellar mass reaches Af*. Without the accretion, contribution to the 
total luminosity is now only by the interior one and the effective temperature is given by r*,cff = (L* /47r cri?^)^/^. 
With higher M^, protostars have larger radii and their birthline shifts toward the upper-right. This dependence on A/* is 
consistent with previous calculations (Palla & Stabler 1992; Norbcrg & Maedcr 2000). 

If a high accretion rate of > 10^"* Mq/yi is indeed achieved in massive star formation, stars as massive as 10 Mq 
have not yet arrived at the ZAMS at the end of the accretion phase. Thus, the presence of such massive PMS stars is 
peculiar to the high scenario. These PMS stars are very luminous with > IO'^Lq, and their effective temperature is 
much lower than that of the ZAMS stars. Alt hough some candidates have bee n reported very close to the ZAMS line, no 
such object has yet been firmly detected (e.g.. lHanson. Howarth &: Contilll997t ). One explanation for the lack of detection 
is the very short KH timescale Ikh of such PMS stars. For example, in the case with = lO~^M0/yr (run MD3), ^kh 
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Fig. 21. — Evolution of the effective temperatures and bolometric luminosities of protostars with accretion rate A/, = 10~^ (solid; run 
MD3), 3 X 10~^ (da shed; MDSxS), and 6 X IQ-^MQ/yr (dot-dashed; MD6x3). The condition required for the protostar in Orion KL region 
IIMorino et al.lll998l) is indicated by the horizontal lines with arrows; Ltot > 4 X 10^ Lq and T^g < 5500K. The mass range satisfying this 
condition is denoted by the shade for the case of M, = 6 X lO-^M0/yr. 



falls below lO^'yr for A/* > lOMg (Fig. [51 bottom). On the other hand, from the case with = 10"^ Mo/yr (run MD5) 
we see that iKH > 10® yr (10^ yr) for low- (intermediate-, respectively) mass protostars (Fig. [T^l bottom). In addition to 
the small number of massive stars, such a short duration in the PMS phase severely limits the possibility of detection. 

As referred to in § [1] some info rmation of du st-enshrouded protostars has been gleaned from the light re-emitted in 
the envelope. On the other hand, iMorino et al.l (fT998f ) reported a unique indirect observation of an embedded massive 
protostar. They observed the infrared light from the Orion BN/KL nebula, which is scattered by dust grains far away 
from the exciting protostar. Direct light from the protostar is blocked by a torus-like structure on the line of sight toward 
us, but likely leaks in the polar direction and is scattered by the nebula. The color temperature of the scattered light 
is evaluated as Teff — 3000 — 5500K from features of the absorption lines. The bolometric luminosity of the protostar is 
Lho\ > 4 X IO^'^Lq. If the scattered light is from the photosphere of the star, this low Tea cannot be explained by a MS star: 
Teff of a ZAMS star with the same luminosity is about 35000K, much higher than the above estimation. iMorino et al.l 
(|1998^ concluded th at a very large stellar radius of > SOOi?© is needed to explain the low effective temperature. Stimulated 
by this result. Nakano et al.l (12000) studied the evolution of protostars with very high accretion rates ^ lO^^A/0/yr using 
a one-zone model. However, they concluded that stellar radius does not exceed 30 Rq and then failed to reproduce the 
observed low effective temperature. In contrast, using more detailed modeling, we have shown that the stellar radius in 
fact reaches as large as ~ lOOM© in the cases of high accretion rate . In Figure [HI we present evolution of photospheric 
temperature and bolometric luminosity for three cases with Af* > 10""^ M^g/yr. This figure shows that the required 
condition is satisfied, for example, with the case of M* > 6 x 10"'^ Mg/yr for Af* ~ 20 — 25 Mq. As a conclusion, the 
observation can be explained if the accretion rate is higher than 4 x 10~^ Mq/yt. The allowed mass range corresponds 
to the duration of the swelling of the protostar by the luminosity wave (see 13. ip . Recall that this swelling is caused 
by very inhomogeneous entropy distribution: the matter near the surface receives a large amount of entropy temporarily. 
Since this is not included in Nakano et al. (2000) 's one-zone model, their model gave a smaller radius than ours. In 
Appendix [Cl we present calibration of parameters for a polytropic one-zone model to include this effect approximately. 
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We have studied evolution of accreting protostars by solving the structure of the central growing stars and the surround- 
ing accreting envelopes simultaneously. Particular attention is paid to cases with high accretion rates of M* > 10~^ Mq /yr, 
which are envisaged in some current scenarios of massive star formation. The protostellar evolution at a high mass accre- 
tion rate of M» = 10"^ Mq/ji ( run MD3 in Table 1) can be summarized as follows; 

The entire evolution is divided into four characteristic phases; (I) adiabatic accretion (M* < 6 Mq), (II) swelling 
(6 Mq<M^< 10 Mq), (III) KH contraction (10 Mq < M^ < 30 Mq), and (IV) main-sequence accretion (M* > 30 Mq) 
phases. 

A main driver of transition between the evolutionary phases above is a decrease in opacity in the stellar interior with 
increasing temperature, and thus protostellar mass. Due to the fast accretion, the matter accreted onto the stellar surface 
is embedded into the interior before radiatively losing the entropy produced at the accretion shock. Early in the evolution, 
the opacity in the stellar interior, which is mainly by free-free absorption, is very high owing to low temperature. As a 
result of the high opacity, the star keeps a large amount of entropy imported by the accreted matter. This is the adiabatic 
accretion phase (I) . With the decrease in opacity, the entropy is transported outward gradually. When the matter near the 
surface temporarily has a large amount of entropy, the star expands as large as ~ 200 Rq. This is the swelling phase (II). 
After that, all parts of the protostar lose energy and the protostar contracts. This is the KH contraction phase (III). With 
contraction, the temperature increases and the hydrogen burning eventually begins. When the nuclear burning provides 
sufficient energy to balance with that lost from the star radiatively, the star stops the KH contraction and reaches the 
zero-age main-sequence (ZAMS) phase at ~ 30 Mq (IV). 

Deuterium burning plays little role in the evolution. Even if the deuterium burning is turned off by hand, the result 
hardly changes. The swelling (II) is caused by outward transfer of entropy within the star, which is observed as a 
propagating luminosity wave. 

In general, at the higher accretion rate, the protostellar radius is larger and then the maximum temperature is lower at 
the same stellar mass, owing to the higher entropy within the protostar. As a result of the lower temperature, the onset 
of the nuclear burning is postponed to the higher protostellar mass. 

For very high accretion rate M* > 3 x 10~^ Mq/jt, in the course of the KH contraction, the radiation pressure onto the 
inner accreting envelope becomes so strong that the flow is retarded before hitting the stellar surface. The reduced ram 
pressure onto the stellar surface causes abrupt expansion of the star. The steady-state accretion is not possible thereafter. 
At this moment, the accretion might be halted, or might continue at a reduced rate or in a sporadic way. In either case, 
the evolution at the critical accretion rate of Md- = 3 x 10~^ Mq/ji gives the upper mass limit of PMS stars at ~ 60 Mq. 
Further growth of the star should be finally limited by the radiation pressure onto the outer dusty envelope. We have 
found this limit with A'/cr gives the maximum mass of MS stars at ^ 250 Mq . 

Such a high accretion rate has also been expected and studied for the first star formation in the universe. The 
evolutions of primordial and solar-metallicity protostars are very similar at the same accretion rate. However, the lower 
opacity of the primordial gas results in the earlier transitions between the evolutionary phases above. For example, with 
M» = 10^'^ Mq /yi, the accreting star enters the KH contraction phase at 10 Mq in the solar metallicity case, while it 
enters at 7Mq in the primordial case. 

Distinguishing features of those massive protostars accreting at high rates include the large stellar radius sometimes 
exceeding 100 Rq, and then low color temperature ~ 6000K. A massive protostar in Orion BN/KL nebula indeed has 
such features and is a possible candidate of such objects. 

The authors thank Francesco Palla, Nanda Kumar, Jonathan Tan, Mark Krumholz, and Shu-ichiro Inutsuka for helpful 
comments and discussions. This stiidy is siipported in part by Research Fellowships of the Japan Society for the Promotion 
of Science for Young Scientists (TH) and by the Grants-in-Aid by the Ministry of Education, Science and Culture of Japan 
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APPENDIX 
A. METHOD OF CALCULATIONS 
A.l. Evolutionary Calculations 
A. 1.1. Protostar 

For protostellar structure, the following four stellar structure equations are solved: 

1 
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(dL_\ ^^_^{ds 



( ds\ _ GM_ fds\ ( L 



\dMj,- 4^r4 \d~p)^\L: V ^' ^^^^ 

with a spatial variable of Lagrangian mass coordinate M . In the above, e is the energy production rate by nuclear fusion, 
s is the entropy per unit mass, Ls is the radiative luminosity with adiabatic temperature gradient, and other quantities 
have ordinary meanings. The coefficient, C in equation (|A4p is unity if L < (i.e., in radiative layers), and given by the 
mixing-length theory if L > Ls (i-C, in convective layers). The momentum equation (IA2p simply states the hydrostatic 
equilibrium. This is valid because the stellar mechanical timescale is generally much shorter than the accretion timescale. 
In energy equation (jA3p . on the other hand, the non-equilibrium term —T{ds/dt)M should be included. This is because 
the KH timescale can be longer than or comparable to the accretion timescale in our calculations. Thus, the thermal 
equilibrium is not generally satisfied in accreting protostars. Nuclear reactions included in e are deuterium burning and 
hydrogen burning via the pp-chain and CN-cycle. The elements are homogeneously distributed in convective layers by 
the convective mixing. 

There are two alternative schemes to integrate equations (jA3[) and (|A4[) . which are the explicit and implicit schemes: 



Explicit Scheme In the explicit scheme, equation (jA3p is considered as a time development equation of the entropy. 
In order to update the entropy at each spatial grid from a time step t — At to t, we use {dL/dM)t evaluated at the 
previous time step t — At and calculate {ds/dt)M by equation (|A3|) . Equation (|A4|) is transformed as, 

Airr^ ^^P\ f ds \ 
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with C = 1 (radiative). We substitute the updated entropy gradient {ds/dM)t to equation (|A5P and calculate the 
luminosity profile L{M,t). This scheme enables the precise integration where the heat transport among mass cells is 
negligible, i.e., {dL/dM)t is much smaller than other terms in equation (|A3p . For example, this quasi-adiabatic situation 
is created in the very opaque interior of the protostar. 

Implicit Scheme In the implicit scheme, on the other hand, time-derivative term —T{ds/dt)M in equation (|A3[) is 
regarded as a source term to calculate the luminosity profile by spatial integration. We simultaneously integrate equations 
(jA3[) and (|A4[) . evaluating {ds/dt)M at every increment of spatial grids. This scheme works effectively where the heat 
transport among mass cells works to redistribute entropy. For example, this is the case near the stellar surface. Devote 
special attention to discretizing the time-derivative {ds/dt)M in equation (jA3p . For example, consider the gas element 
M^,At measured from the surface in a stellar model. This element accretes to the star from a time step t — At to t, 
and is still in the accreting envelope at the previous time step t ~ At. The accreting matter plunges to the star through 
the accretion shock front, which is treated as a mathematical discontinuity in our formulation. Describing the derivative 
across the discontinuity is a problem. We avoid this difficulty adopting a coordinate conversion (SSTSOb), 

s) ^u.m^ m 



where m is the inverse mass coordinate m = Af* — M. The term {ds/dt)m can be obtained only with the entropy 
distribution within the shock front. This treatment is fairly valid, because the structure of the settling layer hardly 
changes over a few time steps, observing from the shock front. 

In our calculations, we properly adopt different schemes in different evolutionary stages. In the early phase of the 
evolution, we use both the explicit and implicit schemes to construct one stellar model (e.g., SSTSOb); the explicit scheme 
is used in the deep interior and is switched to the implicit scheme near the stellar surface. First, we solve the interior with 
the explicit scheme; we integrate equations (jAl[) and (jA2p outward with a guess of central pressure Pc- When the explicitly 
solved region is radiative, we calculate L(M,t) with equation (jA3p . Once active deuterium burning begins, however, the 
explicit update of entropy sometimes makes negative entropy gradients {ds/dM)t < 0. Such a layer should be regarded 



28 



T. Hosokawa & K. Omukai 



as a convective layer, and equation (jASp is not valid any more. We temporarily change the integration method in such a 
layer. We correct the entropy distribution as {ds/dM)t = there after the explicit update, following the "equal-area rule" 
developed by SSTSOb. Once the entropy profile is corrected, we calculate {ds/dt)M and obtain the luminosity distribution 
by integrating equation (jA3p . As the integration approaches the surface, the term {dL/dM)t in equation (|A3p gradually 
grows. We switch to the implicit method where {dL/dM)t attains 0.5AI^,T{ds/dm)t- We integrate all equations (jAip - 
l|A4p with a guess of the entropy at the switching point Sgw When the integration reaches the stellar surface, we check 
the differences from the required boundary conditions (see i; IA.1.2l below). We improve the guessed values of Pc and Ssw, 
and iterate this procedure until the model satisfies the outer boundary conditions. 

The above method successfully works in the early phase of evolution. However, only in the later stage of evolution, 
where the entropy redistribution among mass cells works even in the deep stellar interior, we need to use another method 
with the implicit scheme. This occurs when a widespread convective layer appears during active deuterium burning, or 
when radiative heat transport becomes efficient due to the opacity decrease in the deep interior (see § [3] for detail) . In the 
fully implicit method, we integrate the equations (|Aip ~ (|A4p outward from the center with a guess of the central pressure 
Pc and entropy Sc- If the integration reaches the stellar surface, we advance to the same fitting process at the surface as 
in the above method. If the integration diverges and fails before reaching the surface, we adopt another shooting method 
to a halfway fitting point. We additionally guess the radius i?* and luminosity at the stellar surface, and integrate 
backward to the fitting point. We adjust boundary values to match the physical quantities at the fitting point, which are 
obtained by the forward and backward integrations. 

A. 1.2. Accreting Gas Envelope and Accretion Shock Front 

The outer boundary condition of the protostar is given by constructing the structure of an accreting gas envelope. The 
boundary layer connecting the accreting envelope and protostar actually has detailed structure. The acreted gas initially 
hits a standing shock front and is heated up by compression there. The shocked gas next enters the relaxation layer behind 
the shock front and is cooled down by radiative loss. In our formulation, this detailed structure is not solved, but treated 
as a mathematical discontinuity (SST80). We define the following two points across the discontinuity: a post-shock point, 
where gas and radiation temperature first become equal behind the relaxation layer, and a post-shock point ahead of the 
shock front. In our calculations, the former is the outer boundary of the protostar, and the latter is the inner boundary 
of the accreting envelope. Physical states at each point are related by a radiative shock jump condition. 

Here, we focus on cases where the post-shock quantities are given in advance by the trial outward integration of the 
protostar (see ij lA.l.ip . When the outward integration fails and the shooting to a halfway fitting point is needed, we only 
slightly modify the same procedure (e.g., see PS91). First, we judge if the accreting envelope is optically thin or thick. 
Temporarily assuming that the the accretion fiow is optically thin and in the free-fall, velocity and mass density at the 
pre-shock point are, 



•^pre 



(A7) 



/^p- = J-^—- (A8) 

Ahead of the shock front, outward radiation comes from both the stellar interior and relaxation layer behind the shock 
front. Therefore, radiation temperature at the pre-shock point is written as, 

where a is the Stefan-Boltzmann constant, and Lace is the accretion luminosity defined by equation ([3]) Using equations 
(|A7p - (jAOp . optical thickness of the accreting envelope is estimated as, 

■''pre — -R*Ppre'^(Pprci Tprc); (^10) 

where k is the Rosseland mean opacity. If Tpio < 1, the accreting envelope is optically thin and initial adoption of the 
thin envelope is verified. The jump condition at the accretion shock front is provided by considering the flow structure 
across the shock front (SST80b). In the optically-thin case, the post-shock temperature is related to the stellar radius 
and luminosity (PS91), 
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This is the outer boundary condition of the protostar. The protostellar models are constructed for the post-shock quantities 
to satisfy equation (jAlip . 

If Tprc > 1, on the other hand, the accreting envelope should actually be optically thick. In this case, we solve structure 
of the flow inside the photosphere, i.e., radiative precursor. First, we determine the locus of the photosphere. Assuming 
the free-fall flow outside the photosphere, optical depth to the photosphere Tph can be written following equations (|A7p 
- (jAlOp by substituting physical quantities at the photosphere. The only unknown quantity is the luminosity at the 
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photosphere iph, which is needed to calculate temperature in equations (|A9|) . We calculate this photospheric luminosity 
using the conservation law of the net energy outflow rate Lf, through the precursor, 

= Lph - ih (^V + ^Uph - ^) , (A12) 

where /iph = /i(ppii,Tph) is the enthalpy of both gas and radiation at the photosphere. Using the known post-shock 
quantities, we can give the numerical value of Lg at the bottom of the accretion flow, 

Le^L,- [h^ost + \u%, - ^) . (A13) 

Assigning this value to Le in equation (|A12p . the photospheric luminosity is written with the infall velocity and thermal 
state at r = i?ph. Consequently, Tph and Tph are given as implicit functions of i?ph- Therefore, we can fix i?ph in an 
iterative fashion by requiring Tph = 1. Once the photospheric radius is fixed, we solve the structure of the precursor by 
integrating the following equations inward from the photosphere. 
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where F is the radiative flux. Using equations (jA12p and (|A14p . F is written as a function of p, u, and T, 

F^-pu(^ + h+lu^- ^] . (A17) 

When the integration reaches the pre-shock point, a jump condition bridges the the pre-shock and post-shock points. An 
isothermal jump condition is applied for the optically-thick radiative shock front. 

This is another outer boundary condition in addition to equation (jAlip . 

A. 2. Initial Models 

We construct initial models by solving their structure. The method of calculation is based on that of evolutionary 
calculations, but with slight modifications. First, we assume the entropy distribution in the stellar interior as an increasing 
function with M , 

s(M)^s,,o + f3—^, (A19) 
mH Mh.,0 

where Sc,o is specific entropy at the stellar center, /cb is Boltzmann constant, mn is atomic mass unit, and (3 (> 0) is 
a free parameter. In the core interior, we integrate equations (|Aip and (jA2p outward beginning with guessed Sc,o and 
central pressure, pc- According to equation (jASp . the luminosity profile there is calculated using an entropy gradient 
obtained from (jA19p . Whereas this procedure is valid in the adiabatic interior, entropy and luminosity profiles should be 
consistently calculated in place of equation (jA19p near the surface (also see § lA.l.ip . The switching point is defined as 
follows. Using equations (jA3p and (jA6p . the energy equation can be written as. 



omitting the nuclear energy production and time-derivative terms. Once the luminosity gradient {dL/dM)t becomes 
comparable to M^,T{ds/dM)t calculated with equation (|A19p . we solve the full equations including equation (|A20p . The 
guessed Sc,o and Pc are adjusted so that the core satisfies outer boundary conditions, as explained in § IA.1.21 We have 
confirmed that some variations of initial models only affect subsequent protostellar evolution in a very early phase. In 
this paper, we adopt /3 = 1 as the entropy slope in equation (|A19p . 

A. 3. Opacity Tables 

For the opacity, we adopt OPAL tables (e.g.. Ilgresias fc Rogersl[T996l) with the composition given byjGrevesse & Noeld 

(1993) for temperature T > 7000 K. For T < 7000 K, we use other opacity tables based on calculations by [A lexander & Fcrg usonI 

(1994) . Contribution from grains is excluded in the adopted composition (jAsplund. Grevesse &: Sauva]|200 5: Cmiha, Hubcn v &: Lan: 
I2OOQ) . In all runs, we adopt X = 0.7, Y = 1—X — Z, and relative abundances of heavv elements following fCrevesse fc Noelsi 

(|1993( ) composition. 
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B. COMPARISON WITH PREVIOUS WORK 
- EVOLUTION WITH THE LOW ACCRETION RATE OF 10^^ ^710 /yr - 

In § 13.21 we have presented the calculated protostellar evolution at the low accretion rate of 10~^ Mq/jt (run MD5). 
Whereas our numerical results reproduce basic features obtained by previous work (e.g., SSTSOa, PS91), there are still 
some differences. In this appendix, we examine causes of these differences by exploring some cases with different deuterium 
abundances and initial models. 

B.l. Dependence on Deuterium Abundance 

The effect of different deuterium abundances on the stellar interior structure in early phases is shown in Figure [B22l for 
cases with [D/H] = 3 x 10^^, 2.5 x 10^^, and 1 x 10^^. The middle panel is the same as Figure [7] (upper), but reproduced 
for comparison. As described in 13.21 our run MD5 exhibits somewhat complex evolution of convective layers in the early 
phase. A convective layer first appears at ~ 0.3 Mq, just after the deuterium burning begins. This convective layer 
disappears at ~ 0.6 Mq. The second convective layer emerges just outside the former outer boundary of the first one 
at ~ 0.7 Mq, and quickly reaches the surface. SSTSOa calculated the same model, but did not observe such a complex 
evolution. In their calculation, the first convective layer survives and eventually incorporates most of the stellar interior. 

This difference comes from the fact that the evolution of convective layers sensitively depends on the initial deuterium 
abundance. For example, in the case with the deuterium abundance [D/H] = 3 x 10~^ (Figure [B22I upper: run MD5-dh3), 
slightly higher than the fiducial value, the convective layer first appears at the same epoch as the fiducial case with [D/H] 
= 2.5 x 10~^, and extends outward. In this case, another convective layer soon emerges at the surface and begins to 
extend inward. These two convective layers finally merge at ~ 0.6 Mq, and render most of the interior convective. At 
M^ ~ 1 A/0, the convective region contains about 97% of the total mass. These features are very similar to those found 
by SSTSOa, who terminated their calculation at Af* = 1 Mq. Our calculation shows that the subsequent evolution in 
this case is almost the same as in our fiducial run MD5. The central radiative core gradually expands as the stellar mass 
increases, and finally occupies most of the interior at ~ 4 Mq. 

This sensitive dependence on the deuterium abundance is explained as follows. The expansion of convective layers is 
caused by increase of the specific entropy within it, which is generated by the nuclear burning at the bottom of the layer. 
For continuing expansion of the convective layer, deuterium needs to be supplied continuously. Fresh deuterium becomes 
available when the convective layer newly incorporates the radiative region. This deuterium immediately spreads over 
the convective layer, and is consumed for further nuclear fusion. For sufficiently high deuterium abundance, acquired 
deuterium is high enough to maintain the active nuclear burning. In this case, the convective layer continuously grows as 
in the case with [D/H] = 3 x 10^^ (run MD5-dh3). On the other hand, for low deuterium abundances, the nuclear burning 
becomes too weak to maintain development of the convective layer. The convective layer ceases to grow and disappears 
as in the case with [D/H] = 2.5 x lO^'^ (run MD5). For even smaller deuterium abundance [D/H] = 1 x 10~^ fFig. IB22"1 
lower; run MD5-dhl), no convective layer appears even after ignition of the deuterium burning until at Af* ^ 0.8 Af0, 
when the surface convective layer finally emerges. In this case, the subsequent mass-radius relation is almost same as in 
the case without deuterium burning (run MD5-noD) although the surface convective layer pesists until M^, ~ 3 A'/0. Our 
fiducial run MD5 is just intermediate between the deuterium-rich (run MD5-dh3) and -poor (MD5-dhl) cases presented 
in Figure [B22I Therefore, a slight difference in input physics affects the internal structure significantly. 

B.2. Effect of Difference in Initial Models 

Evolution of accreting intermediate- mass (1 Af0 < A/* < S Mq) protostars have been studied in some previous work 
(e.g., PS91). Our results agree well with the previous work, but still show some differences even with the same accretion 
rate, which can be ascribed to different initial models. Whereas our initial model is a tiny radiative protostar with ^ IMq 
(see ij IA.2[) . the initial model of PS91, for example, is a 1 Mq fuUy-convective protostar. We confirmed that our code 
reproduces almost the same evolution as that of PS91, if we begin the calculation with the same initial model. 

The initial model of PS91 is constructed in a manner similar to that described in § IA.2I The homo geneous e ntrop y 
distribution is adopted in place of equation (|A19I) . and its value is taken from the semi-analytic model of IStahle"rl ()19S8[ ). 
We have reconstructed the model of iStahlen (|1988f) . and obtained Sc,o = —4.12 /cb/'tih for a 1 Mq star with our adopted 
opacity tables. As in ii lA.2[ we separately solve the core interior and the surface super-adiabatic layer. First, we integrate 
the core interior outward from the center with a guessed value of the central pressure. We also guess the luminosity at 
the bottom of the surface super-adiabatic layer in advance, and record the entropy gradient ds/dM using equation (jA4p . 
Once the recorded ds/dM attains a small finite value, we switch the scheme and integrate all equations (jAip - l|A4p . The 
guessed quantities are adjusted for the core to satisfy the outer boundary conditions. The constructed model is a fully 
convective star, whose mass and radius are 1 Mq and 4.2 Rq. 

We calculate the subsequent evolution beginning with this initial model (run MD5-ps91). Figure [B23l shows the evolution 
of such a protostar. Our calculation reproduces the basic features of a calculation done by PS91, which adopted the same 
accretion rate and accretion-shock outer boundary. In an early phase, the star remains fully convective and deuterium 
burning occurs around the stellar center. Deuterium concentration significantly falls in this phase. This is because the 
central temperature c ontinuously increases, and total burning rate Ld slightly exceeds the steady burning rate LD,st (e.g., 
iPalla fc Stahia1ll993D . These features are somewhat different from our fiducial run MD5, where a radiative core persists 
throughout the evolution and gradually expands with the increase in A/* by accretion (see Fig. [7]). These differences 
come from different initial models. Figure [7] also shows that mass-averaged deuterium concentration /d,av does not fall 
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much below 0.01 in our run MD5. This is because the deuterium burning layer gradually moves to the outer layer, where 
temperature is always around 10® K and Ld ^ -^D.st- At M^, — 2 Mq, for example, /d,av ~ 0.04 in run MD5, but only 
/d,av ~ 10~^ in run MD5-ps91. Therefore, our calculation predicts that intermediate-mass protostars have a larger amount 
of deuterium than assumed by PS91. Figure [B23l shows that the fully convective phase ends at ~ 2.6 Mq, when a thin 
radiative layer, i.e., radiative barrier, suddenly appears within the star. The radiative barrier blocks inward convective 
transport of the accreted deuterium. Consequently, the inner region, where deuterium has run out, immediately returns 
to being radiative. After that, the convective layer remains only above the radiative core, and deuterium burning occurs 
at the bottom of the convective layer. As in our fiducial run MD5, the protostar swells up at ~ 3.4 AIq. PS91 have 




logM,(Me) 

Fig. B22. — The interior evolution in early phases for different deuterium abundances with the accretion rate A/» = 10~^ Mq/jt. The 
deuterium abundances are [D/H] = 3 X 10"^ (top; run MD5-dh3), 2.5 X 10~^ (middle; MD5), and 1 X 10"^ (bottom; MD5-dhl). Properties of 
the interior structure are presented in the same manner as in Fig. |2] In each panel, the thick solid line represents the position of the accretion 
shock front. The convective regions are shown by the gray-shaded area. In the top (bottom) panel, the position of the accretion shock front 
for [D/H] = 2.5 X 10~^ ([D/H]=0, run MD5-noD) is plotted with the dot-dashed curve for comparison. 
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logM.CMe) 

Fig. B23. — Same as Fig. [7] but for the different initial model calculated following PS91 (run MD5-ps91). In the upper panel, the filled 
circles denote the protostellar radius by PS91 with the same parameter as our fiducial case (MD5), which is shown by the dashed curve. 
PS91 started calculation at M.^o = 1 A^q, while the initial mass in our calculation is much smaller M,^o = 0.01 Mq. The dot-dashed curve 
also presents the protostellar radius in run MD5-ps91, where the deuterium burning is turned off at Af* = 2.7 AIq, when a radiative barrier 
emerges within the star (see text). 



concluded that this swelling is caused by the shell-burning of deuterium. To see the role of the deuterium shell-burning 
in the swelling, we also calculate evolution with cessation of the deuterium burning after the appearance of the radiative 
barrier. We find that the swelling occurs similarly, although it is slightly delayed, even without deuterium burning. This 
is caused by the outward transport of embedded entropy, which is observed as propagation of a luminosity wave (see 
ij 13.11 and l3.2p . Although PS91 have attributed the swelling only to the shell-burning of deuterium, we conclude that the 
main driver of the swelling is always the propagation of the luminosity wave. After the turn-around of the radius, the 
star enters the subsequent KH contraction phase. The calculated mass-radius relation thereafter gradually converges with 
that in our fiducial run MD5. The epoch of the Hydrogen burning, ~ 7 Mq, also agrees with that in the fuducial run. 

C. CALIBRATION OF ONE-ZONE MODEL 

While detailed numerical calculations as in this paper are a direct method to study the structure of accreting protostars, 
one-zone modeling with the polytropic equation of state P — Kp^~^^^^ is also a useful way. In the one-zone models of 
protostellar evolution, which were originally developed by Nakano, Hasegawa, & Norman (1995), an energy budget of 
the protostar is considered. The total energy of the star E and its time-derivative dE/dt are written as functions of 
the stellar mass and radius i?* respectively. The mass-radius relation is obtained by substituting the equation E 
into dE/dt. iNakano et a l. (2000) applied such a model to study protostellar evolution at the very high accretion rate 
of M* 10"^ Mq/jt. iMcKee fc Tanl (|2003t ) improved the model by cahbrating it with numerical calculations at the 
accretion rates of < lO"'^ Afg/yr (e.g., PS91). iTan fc McKe3 (|2004D also showed that the one-zon e models can reproduc e 
the mass-radius relations of primordial protostars at the high accretion rates lO^'^ Mq/jt fe.g.. lOmukai fc Pallall2003[ ). 
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Here, we pre s ent on e-zone models reproducing our numerical results fFig. lC24|l . These models are constructed following 
iMcKee fc Tal (l200l with some improvements. Our numerical calculations have shown that protostars are initially 
radiative for any accretion rate. With an accretion rate less than 10~^ Mq/jt, the convective layers gradually extend 
as deuterium burning becomes significant. Since the pol ytropic m o dels w ith index N = 1.5 and = 3 approximate 
a fully convective and radiative star respectively (e.g.. ICox fc Giulilll968[ ). we adopt intermediate values: N = 2.5 for 
< 10~^ Mq/jt and shghtly larger values for higher accretion rates, = 2.5 + 0.25 log(M,/10~^ Mo/yr). Initial 
masses and radii are arbitrarily taken as M<,^o = 0.1 Mq and R^^ = 2.5 i?o(Af,/10^^ Mo/yr)^/^. The steady deuterium 
burning is assumed to begin when the central temperature reaches 1.5 x 10^ K. After that, we adopt N — 1.75, which is 
slightly larger than the fully convective value N = 1.5, to resemble our numerical result, where a radiative core gradually 
expands with the increase of the stellar mass (see § 13.21) . Subsequently, the protostar swells up with the propagation of 
the luminosity wave. We include this effect into the model by artificially increasing the stellar radius by a factor of 3 to fit 
the calculated mass-radius relations. We assume that this occurs when the ratio iKu/^acc reaches 1.75. Other adjustments 
of models are the same as those of iMcKee fc TanI ()2003f) . The presented one-zone models fit our mass-radius relations 
within 30% deviation except in very early phases for high accretion rates. 



